The Annals of Applied Probability 
2009, Vol. 19, No. 4, 1319-1346 
DOI: 10.1214/08-AAP571 

© Institute of Mathematical Statisties, 2009 



ON THE INVERSE FIRST-PASSAGE-TIME PROBLEM 
FOR A WIENER PROCESS^'^ 

By Cristina Zucca and Laura Sacerdote 

University of Torino 

The inverse first-passage problem for a Wiener process {Wt)t>o 
seeks to determine a function b : R+ — ^ R such that 



1. Introduction. Many phenomena can be modeled as first passage time 
of suitable Markov processes through constant or time varying boundaries. 
The first-passage problem has a long history and a large number of appli- 
cations that range from finance to engineering and biology [see, e.g., Ric- 
ciardi and Sato (1990) for some references]. Yet explicit solutions to the 
first-passage problem are known only in a limited number of special cases 
(including linear or quadratic boundaries). 

In modeling generally one describes the dynamics of the involved vari- 
ables via a suitable stochastic process {Xt,t > 0} constrained by an assigned 
boundary b: {0,oo) M satisfying 6(0+) > and investigates distribution 
features of the first passage time (FPT) 
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T^mf{t>0\Wt>b{t)} 



has a given law. In this paper two methods for approximating the 
unknown function b are presented. The errors of the two methods are 
studied. A set of examples illustrates the methods. Possible applica- 
tions are enlighted. 




n = mi{t>0\Wt>b{t)} 
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of Xt over h. This is the direct FPT problem. However, there are also in- 
stances when the underlying stochastic process is assigned, one knows or 
estimates the FPT distribution Fi, and wishes to determine the correspond- 
ing boundary shape. This is the inverse first passage time (IFPT) problem. 
Different applications of the IFPT problem can be listed. In a reliability con- 
text one could compare performances of alternative objects, characterized 
by their lifetime distribution comparing the shapes of the boundaries. In this 
case one should assume the same underlying stochastic process to describe 
the evolution toward the crash of the object and to identify the crash time 
as an observed FPT of the process through an unknown boundary. Other 
interesting applications can arise when the resulting boundary exhibits a 
periodic behavior. This fact can suggest particular features of the modeled 
phenomenon in a finance or in a neurobiological context. Despite its impor- 
tance in applications, the literature focuses only on specific problems [cf. 
Sacerdote and Zucca (2003a, 2003b) and Sacerdote, Villa and Zucca (2006)] 
while there is a lack of mathematical results. Here we try to cover this gap 
focusing on the inverse first passage time problem for the Wiener process. 

The IFPT problem was brought to our attention by Professor Goran 
Peskir, who presented us its first formulation by A. Shiryaev in 1976 (during 
a Banach center meeting). The original Shiryaev question considered the 
case when Fi, is an exponential distribution. There exist two early papers 
dealing with the existence problem written by Dudley and Gutmann (1977) 
and Anulova (1980). These papers provide the existence of some stopping 
times for given F^; however, these stopping times are not of the form (1.1) 
for some function b. Furthermore, in Capocelli and Ricciardi (1972) the 
properties characterizing the FPT distribution in the case of the constant 
boundary are discussed. 

The problem of the existence and uniqueness of b is still open. Through 
the paper we do not enter into the question of the existence of b but assume 
that such a boundary exists, is unique and sufficiently regular (continuous 
or C^). 

The main aim of the present paper is to propose two algorithms for ap- 
proximating the unknown function b when F^, is given. After some prelimi- 
naries presented in Section 2, in Section 3 we introduce a first approach to 
the IFPT problem. It is based on the idea to approximate 6 by a piecewise- 
linear boundary c for which the law of Tc can be computed. It leads to a 
Monte Carlo method for determining c. The second approach (Section 5) is 
based on the classic idea due to Volterra (around 1896) on how to approx- 
imate an integral equation (of the first kind) by a system of n equations 
in n unknowns. Since the integral equation for b is nonlinear, the resulting 
system is also nonlinear. We propose a numerical method for approximating 
b at finitely many points. 
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Both methods produce an approximation of the boundary value. Hence it 
is necessary to evaluate the respective errors, that is, the difference between 
the exact unknown value and the computed approximation. We limit the 
study of the error at the discretization time knots obtained with constant 
discretization step h. The error of the first method is discussed in Section 4. 
It is due to multiple causes but it is dominated by the error due to the 
substitution of the continuous boundary b with a piecewise-linear boundary 
c and it is 0{h'^) . Furthermore, due to the use of the Monte Carlo method to 
evaluate some involved integrals, the resulting boundary is estimated with 
a fixed confidence level a. The error of the second approach is evaluated in 
Section 6 and it is 0(h) or 0{h?) depending upon the numerical method used 
to evaluate the involved integral. In Section 7 we illustrate the goodness of 
the proposed methods by means of a set of examples. Finally, in Section 8 we 
consider the boundary corresponding to the exponential FPT distribution 
and we give a numerical answer to the 1976 Shiryaev question. 

2. First passage time: the direct and inverse problem. Given a stan- 
dard Wiener process W = {Wt)t>o started at zero, and a sufficiently regular 
function 6: (0, oo) M satisfying 6(0+) > 0, denote 

(2.1) Tb = inf{t >0|Wt >6(t)} 
the FPT of W over 6, and set 

(2.2) f^it) = ±W(r,<t) 

to denote its density function for t > 0. In the sequel we also need to consider 
the Wiener process starting in xq at time to; in this case we write fb{t\xQ,to) 
in spite of fb{t)- Throughout we assume that all regularity assumptions 
ensuring the existence of the objects introduced and properties imposed are 
fulfilled. 

By pa^b{t-,x\s,y), we denote the transition probability density function of 
W ait constrained by the absorbing boundary b over [s, t] given that Wg = y, 
that is, 

(2.3) pa,b{t.x\s,y) = ^nWt <x,n> t\Ws = y) 

for X < b{t) with t > s > and y < b{s) given and fixed. 

Under the hypothesis that the stochastic process (in this case the Wiener 
process) is assigned, the direct first-passage time problem seeks to determine 
Ff, when b is given. The inverse first-passage time problem seeks to deter- 
mine b when Ffj is given. It is interesting to note that the IFPT problem 
looks for the boundary b{t) that is a deterministic time dependent function. 
Furthermore, as proved by Strassen (1967), if the boundary b G C^(M'''), 
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then the probabihty distribution of the FPT is absolutely continuous with 
continuous density. 

In the literature there exist some equations that relate quantities (2.2) 
and (2.3), allowing, in a limited number of cases, to solve the direct FPT 
problem [cf. Ricciardi and Sato (1990)]. In Section 2.1 we list some exist- 
ing closed form results for the Wiener process. These results will be used 
in Sections 3 and 7 to numerically validate the reliability of the two algo- 
rithms proposed. In absence of the analytical solution, the direct problem 
can be solved numerically, making use of one of the algorithms proposed in 
the literature [Buonocore, Nobile and Ricciardi (1987) and Zucca (2002)]. In 
Section 7 we use the algorithm introduced in Buonocore, Nobile and Riccia- 
rdi (1987) to estimate the FPT p.d.f. for a set of assigned boundaries. These 
numerical evaluations make it possible to enlarge the test set for the valida- 
tion of the IFPT algorithms proposed in this paper. Finally, we briefly recall 
in Section 2.2 some results that will be useful later on [cf. Peskir (2002)], 
since one of the two algorithms proposed for the IFPT method requires the 
knowledge of 6(0). 

In the sequel we will assume that to = ^^o = when this gives no loss of 
generality. 



2.1. Known boundaries. The first-passage time density fi, for a Wiener 
process is known explicitly in a few cases. The following three will be of 
interest in the sequel: 

1. Linear boundary. If the boundary is given by 
(2.4) c{t)=a + l3t 

with a > xq and /3 G M, then [cf. Doob (1949), page 397, and Malmquist 
(1954), page 526] 

" - ^0 ^-(a+/3(t-fo)-xo)V(2(f-fo)) 

v/27r(t-to)^ 
a-XQ i'a + P{t - to) - xo\ 

V^W^'^y vit-to) ) 

for t > 0, where ip{x) = (l/\/27r)e~^'^/^ is the standard normal probability 
density function. Note that (2.5) is known as the Bachelier-Levy formula. 

2. Daniels boundary. If the boundary is given by 

(2.6) d{t) = f - ^ log + ^/^T^^) , 



fc{t\to,Xo) 



(2.5) 
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where a > 0, /3 > and 7 > -/^V^, then [cf. Daniels (1969)] 
(2.7) 



d(tf/{2t) _ P-(d{t)-af/(2t) 



a 



1 (,Jd{t)\ (3 (d{t)- 
for t > 0. 

3. Piecewise-linear boundary. If the boundary is given by 

(2.8) c{t) = ai + (3it for t G and i > 1, 

where ti = tQ + ih, h> and Oj, /Jj G M with Iq > 0. Setting Oj+i = + /3jti, 
we get that t ^ c{t) is continuous on [to, 00). Let us denote by q = c{ti) the 
knots of c for i > 0. 

The transition density function of in xi, X2, . • . , x„ at ti,t2, ■ ■ ■ ,tn con- 
strained by the absorbing piecewise-hnear boundary c over [Iq , t„] given that 
l^to = 2;o < is [recah (2.3) above] 

Pa,c {ti,Xi;t2,X2; ■ ■ .■,tn,Xn\to,Xo) 

n 

— J_ J_ Pa,citi , Xi\ti—i, a^j—i ) 



(2.9) 



i=l 
n 



e 



-2{ci-Xi){ci-i-Xi-i)/{ti-ti-i)\ 



exp 



(Xj Xj_l) 



for (xi, X2, . . . , Xn) G M" with xi < ci for 1 < i < n and xq < cq where to < 
ti <t2 < ■ ■ ■ <tn are given and fixed. This imphes 

nWt, eCi,...,Wt„&Cn,T,> tn\Wto = Xo) 

(2.10) 

Pa,c(ii,2;i;...;t„,Xn|to,xo)(ixi ■ ■ ■ dxr 



Ci J a 

for any Borel set Ci C (—00, Cj] with 1 < i < n. 

The identity (2.9) is a weh-known fact in the case n = 1 [cf. Doob (1949), 
equation (4.2) and Section 5, or Durbin (1971), Lemma 1]. The case of 
general n > 2 [cf. Wang and Potzelberger (1997)] follows readily by induction 
arguments using that W has stationary independent increments. 

Since there are only a few cases where the FPT density function is known 
in closed form, there has been a big effort in the past to find numerical 
approximation for it. One of the most used algorithms has been presented 
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in Buonocore, Nobile and Ricciardi (1987). This method solves a Volterra 
integral equation of the second kind derived from the Fortet equation but 
characterized by a nonsingular kernel. In Section 7 we apply this method 
to determine numerically the FPT density function in the case of a time 
periodic oscihating boundary. 

2.2. Limits at zero. One of the key issues in the numerical treatment of 
the inverse first-passage problem is to know 6(0+) in terms of /b(0+) and 
vice versa. In this section we display some known results on this relation. 
We will make use of these facts in Sections 7 and 8 below. 

In the notation of Section 2, let W = {Wt)t>o be a standard Wiener pro- 
cess started at zero, let b:{0,oo) ^ M be a continuous function satisfying 
6(0+) > 0, and let us assume that the first-passage time Tb from (2.1) has 
a continuous density function having a limit /fe(0+) in [0, oo]. Then the 
following facts are known to hold [cf. Peskir (2002)]. 

If there are e > and 5 > such that 



(2.11) b{t)>^i2 + e)tlog{l/t) 

for all t £ (0, 6), then fbiO+) = 0. If there is 5 > such that 



(2.12) b{t)<^2tlog{l/t) 

for all t £ {0,5), then /b(0+) = +oo. Moreover, if we define a boundary g by 
setting 



(2.13) g{t) = ^2t log(l/t) + t log Iog(l/t) + ct 

for t G (0, 6c) with c G M given and fixed, then the limit fg{0+) exists and is 
given by 

p-c/2 

(2.14) fg{0+) = 



Observe that g from (2.13) locally at zero lies between the two functions 
appearing on the l.h.s. of (2.11) and (2.12) respectively (where e > may 
be as small as one likes). 



3. First approach: a piecewise linear approximation via Monte Carlo 
simulation. In this section we face the IFPT problem by looking for a 
piecewise-linear approximation of the unknown boundary b with an approach 
that can be considered analogous to the method used by Durbin (1971) for 
the direct FPT problem. Thus, we determine the piecewise-linear boundary 
c (2.8) that approximates the exact boundary b corresponding to the given 
first-passage density /;,. 
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First we consider the inverse problem when /fe(0+) = and we assume 
that 6(0+) > is known. Then, using the notation of Section 2 (with to = 0), 
we set Co = b{0+) and, thus, ai = 5(0+) as well. When 6(0+) is unknown 
we can guess its value, but in this case the error at the initial knot can 
dominate the final one. Since the resulting algorithm uses the Monte Carlo 
approach to estimate involved integrals, henceforth we call it the PLMC 
(Piecewise Linear Monte Carlo) method. In the following we consider a time 
discretization ti = to + ih, i = 1, 2, . . . , where /i is a positive constant. The 
choice h constant is made to simplify our notation, but the method can be 
easily extended to a nonconstant h. 

The driving idea of our algorithmic approach is to equate the probabil- 
ity of the FPT of W through c in {ti-i,ti\ and the analogous probability 
determined by the given first-passage density /{, in the same interval. The 
resulting equations allow to determine the coefficients Oj and Pi in (2.8) suc- 
cessively on the intervals {ti-i,ti] for i > 1. The algorithm can be divided in 
two successive steps. 



Step 1. We determine the value of Pi in (2.8) such that the probability 
of the first-passage of W through c in (0,ti] equals the probability of the 
first-passage of W through b in {0,ti], that is, we look for the value of P 
such that 

(3.1) r ut)dt= r h{t)dt, 



10 JO 

where fc{t) is given by (2.5). 

Step 2. Given ai, . . . , and Pi, . . . ,Pn with n > 1, we set On+i = + 
Pntn and we determine the value of Pn+i such that the probability of the 
first-passage time of W through c in (t„,t„+i] equals the probability of the 
first-passage time of W through b in (t„, tn+i]? that is, we look for the value 
of P such that 

/ / ■■■/ ( fc{t\tn,Xn)Y[pa,c{ti,Xi\ti-i,Xi-iU dxi- ■ ■ dXndt 

(3.2) 

/•l.n + 1 

h{t)dt, 



1=1 



where fc{t\tn, Xn) is the density function of the first-passage time of W over 
On+i + Pn+it for t > tn given that Wt^ = x„, that is, it is given by (2.5) 
with a := c„ and P = Pn+i, while pa,ciii,Xi;ti-i,Xi-i) is given by (2.9) for 
1 <i <n (where to = xq = 0). 
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Remark 3.1. The product appearing on the l.h.s. of (3.2) involves the 
known functions (2.9) depending upon ai, . . . , and (3i,. . . ,(3n determined 
in the preceding step. The unknown On+i and Pn+i appear in (3.2) only 
within the function fc{t\tn,Xn). Furthermore, a^+i can (by continuity of c 
at tn) be expressed in terms of Pn+i and the known a„ and (3n as follows: 

(3.3) On+l = an + iPn - Pn+l)tn 

SO that we only need to determine Pn+i from the equation (3.2). 



Let us now detail the two steps of the algorithm. In order to solve the 
equation (3.1) and (3.2) for the unknown (3i and /?„ respectively, it is nec- 
essary to compute the integrals involved. 

Discretizing the integral of the l.h.s. in (3.1) by a rectangular method 
[cf. Atkinson (1989)], we obtain a nonlinear function of the unknown (3i 
while the r.h.s. of (3.1) is easily computable by standard means. Various 
approximate methods can then be used to solve the resulting equation. Here 
we use the middle point method. 

In the successive steps, mainly when n becomes large, the multiple inte- 
grals of the l.h.s. of (3.2) cannot be handled by standard numerical meth- 
ods. To handle the problem, we adapt the Monte Carlo method proposed 
by Wang and Potzelberger (1997) to our case. 

For this, note that, using expression (2.5) for fc{t\tn,Xn) and (3.3), we get 

tn + l 



rtn + 1 

H{(3n+i;Xn) := fc{t\tn,Xn)dt 

J tn 



Cm Xm / Cn Xn ~\~ (^n+lit ^rt/ i 7, 

V 1= dt 



tn (t-tn)'^/^ \ 

Pn+l{tn+l — tn) + (Cn — Xn) 



1_$ 



■\/tn+l t^ 

(34) I ^-20„+l(cu-Xn)q.f f3n+l{tn+l - tn) - jCn - Xn) \ 

\ \/tn+l — tn ) 

'/5n+l(*n+l — tn) + ttn + f^ntr. 



1-$, 

_j_ g — 2/3„+i(Q:„+/3„f„— x„) 

/3n+l(^n+l ^n) Q^n Pntn ~l~ Xn 



X $ 



■\/tn+l ~ tn 



for each n > 1. Here $(x) = ^{z) dz is the integral of the standard nor- 
mal density. It follows by (2.9) and (3.4) that the equation (3.2) can be 
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rewritten as follows: 

ElH{/3n+i]Xn)f[liXi < c.)e-2(c-^»)(<=»-i-^i-i)/(*-*i-i) 

(3.5) 

= / Mt)dt, 

where Xi ~ N{0,ti — tj-i), i = 1, . . . ,n, are independent random variables 
(and Xq = 0). a Monte Carlo method can now be used to estimate the l.h.s. 
of (3.5). 

In order to find an approximation of /3n+i) we can use, as in the Step 1 
above, the iterative middle point procedure. Using then (3.3), one obtains 
Un+i and thus determines c on {tn,tn+i]- 

Remark 3.2. For each n > 1 given and fixed, the equation (3.5) defines 
a nonlinear function of the unknown parameter (3n+i depending on it only 
through H from (3.4). Since the l.h.s. of (3.5) is monotone in Pn+i, equation 
(3.5) admits a unique solution Pn+i G ^■ 

Remark 3.3. Finally, we modify the algorithm to include the case 
/b(0+) > 0. In this case the equation (3.1) cannot be used on {0,ti] since 
the piecewise-linear boundary c must satisfy c(0+) > 0, while the condition 
fb{0+) > implies c(0+) = (cf. Section 2.2). However, we can use (2.14) to 
estimate the boundary in a neighborhood of zero. The use of this expression 
on the first discretization interval allows then the iteration of the previous 
algorithm on the successive intervals if we disregard the possible crossing 
happened in the neighborhood of zero. 

4. PLMC method error estimation. The errors involved in the approx- 
imation of the boundary values via the PLMC method can be classified in 
three types according to the causes generating them. The first one is due to 
the piecewise linear approximation of the boundary, the second one is due 
to the Monte Carlo approximation of the integrals on the l.h.s. of (3.2) and 
the third one is due to the root method used to compute the zeros of (3.2). 
The last error can be disregarded since it depends on a tolerance that can 
be fixed in advance to make it negligible with respect to the other errors. In 
the following subsections we successively focus the study on the other two 
errors involved in the method. We explicitly underline how these two errors 
have different mathematical natures. Indeed, the first one, as the root error, 
is purely numerical, being determined by our approximation of the unknown 
function via a piecewise linear one. On the contrary, the second error has a 
statistical nature since it is related with the evaluation of an integral via a 
Monte Carlo method. 
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4.1. Error due to the piecewise linear approximation. The PLMC method 
estimates the boundary values through a piecewise hnear approximation, but 

we hmit the study of the error at the knots tn, n = 1, 2, Furthermore, at 

this step we assume the absence of other causes of error and we define error 
of the method at the nth knot the distance between the boundary b and its 

approximated value c on the nth knot |ej^^^*-^| = \b{tn) — c{tn)\, n = 1,2, 

In order to gain an estimate of such error, we first consider in Lemma 4.1 the 
error due to the boundary linearization on a single step. In doing this anal- 
ysis we hypothesize to know the true value of the boundary on the previous 
intervals. A second step considered in Lemma 4.2 studies the propagation 
of the error, that is, we admit an error 6 in the estimation of the boundary 
value at node n — 1 and we evaluate its consequences on the next node n. 
Finally, in Theorem 4.3 we determine the global error of the PLMC Method. 

Lemma 4.1. Let {Wt)t>o be a Wiener process bounded by a monotone 
concave (or convex) boundary 6 G C^([0, oo)). If we approximate in (0, t„] the 
boundary b{t) with the boundary Cn{t), for n = 1,2, . . . , defined by 



(4.1) Cn{t) 



{h{t), ^ tG(0,t„_i], 

1 6(t„_l) +/3„(t - t„_l), t e [tn-ltn], 



the resulting error at the knot as /i — > 0, is 

\b{tn)-c{tn)\^0{h^). 

Proof. We limit ourself to the case of a concave boundary since the 
convex case can be dealt in a similar way. Moreover, we split the proof in 
two parts, first taking into account the first discretization step and later a 
generic step. 

Step 1. Let us consider the following three straight lines (cf. Figure 1) 
on (0,/i]: 

1. Bi{t): the tangent to b{t) in t = 0, 

y = Qi + 6'(0)t. 

2. c(t): the linear boundary determined via the PLMC method in the first 
discretization interval (0, /i], 

y = ai+[3it. 

3. B2{t): the secant through (0,ai) and {h,b{h)), 

b(h) — a\ 
y = ai + ^^ — 
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Due to the concavity and monotonicity hypotheses on the boundary b{t), 
it holds for all t £ (0, h] 

B2{t)<b{t)<Bi{t), 
that implies [cf. Sacerdote and Smith (2004)] 

P(rB, G (0, h]) < P(Te G (0, h]) = F{n G (0, h]) < F{tb, G (0, h]). 
Hence, we get 

B2{t)<cit)<Biit), 
that implies the ordering of the slopes 

b{h) — ai 



(4.2) 



h 



<Pi<b'{0). 



Note that /3i is a function of h: Pi{h). When h^O, due to the hypothesis 
ai = 6(0), the inequality (4.2) implies 

(4.3) /3i=/3i(/i)— .6'(0). 

On the first knot the distance between b{t) and the linear boundary of the 
PLMC method results in 

5"(e) 



\{ai + Pi{h)h)-bih)\ 



(6(0) + Pi{h)h) - (^6(0) + b'{0)h + ^/i') 
{P,{h)-b'mh-^lh' 

no,. 



<\{P,{h)-b'{0))h\ + 



2! 
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<\li'i{vW\ + 



2! 

= 0{h^) + 0{h^) = 0{h^), 

where we made use of (4.3), of the mean value theorem, of the triangular 
inequality and of the McLaurin expansion of b{t) and /3i(t). Here £ (0, /i] 
and r] £ (0, h]. 

Step n. The boundary (4.1) coincides with the boundary b{t) on (0, tri_i) 
and it is determined solving 

(4.4) F{n G [t„-i,tn]) =P(r-G [t„-i,t„]) 

on [tn-i,tn]. Equation (4.4) implies 

nn G [tn~l,tn]\Xt < b{t),t< t„_l}P{Xt < b{t),t< tn-l} 

(4.5) 

= P{t^G [tn^l,tn]\Xt <Cn{t),t< t„_l}P{Xt < C„(t),t < 

and, using (4.1), one has 

r{n€[tn^l,tn]\Xt<b{t),t<tn^l) 

(4.6) 

= P(T,(,_,)+ft^(t_t„_,) G [tn-l,tn]|Xt <6(t),t <t„_i). 

The proof of Step 1 can now be repeated to the conditional probabilities in 
(4.6) in order to complete the proof. □ 

Let us denote On the estimate of a„ and let fin{(in) be the corresponding 
estimate of the slope f3n{an) obtained with the PLMC method on the interval 

Lemma 4.2. An error \5\ = |a„ — a„| propagates on the estimate /?„(a„).' 

\APn\ = \Pnian)-MSn)\-0{6). 

Proof. We separate the proof in two parts; the first one concerning 
the case of the first interval and the second one concerning the generic nth 
interval. 

Step 1. Let us first prove that Pi = /3i(ai) is a continuous and differ- 
entiable function. When ai = 6(0) define 

(1 _ g-2ai/ti(ai+/3itl-^))^^g-^V(2tl) dx - h, 
-oo \'2TTtl 

where ki is the r.h.s. of (3.1). Recognizing the first two terms as the l.h.s. of 
(3.1), equation (3.1) becomes -Fi(ai,/3i) = 0. This is an implicit equation in 
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f3i and ai that admits a continuous and difFerentiable solution Pi = Pi{t) if 
Dini's theorem holds. The hypothesis of Dini's theorem is verified since 



dFi{ai,Pi 



-2ai/ti{ai+l3iti-x)-x^/(2ti) 



d(3i J-oo \/27rt7 

is nonzero for all ai,/3i, having for the hypothesis ai ^ 0. Hence, we get 
that /?! = is a continuous and differentiable function and the error 

connected with the use of Si in spite of ai is propagated on (3i . Hence, /3i 
becomes /3i = /?i + A/3i with 

A/3i=/3i(ai + 5)-/3i(ai)~0(5) 

since /3i(ai) is continuous and differentiable. 

Step n. For n = 2, 3, . . . , we proceed in analogy with step one showing 
that (3n = (3n{cen) IS a continuous and differentiable function. Let us consider 
the approximated stepwise linear boundary c{t) in the time interval [0,t„]. 
Let us define 



Jo Jo 



n-1 

dx\ ■ ■ ■ dXn—l fa{xi, ti\xi--l, ti—i] 
i=0 



X 



1 



(1 - exp{(-2(c„_i - Xn-i) 

X [an + Pntn - Xn)) 
/{tn-tn-l)]) 



X exp 



\Xn Xn—l) 
{2{tn-tn-l)) 



{j27r{tn-tn-l))dXn 



h 



SO that equation (3.2) reads Fn{an, Pn) = 0. This is an implicit equation in 
j3n and a„ that admits a continuous and differentiable solution /3„ = f3n{t) 
if Dini's theorem holds. 
Note that 

dFn{an,f3n) 



df3n 



Jo 



n-1 



fa{xijii\xi—ljt 



i-1) 



1=0 



X / 2tn{Cn—l — Xfi—l) 
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X g-2(c„_i-x„_i)(a„+/3„t„-x„)/(t„-t„_i) 
X g-(a:n-a;„_i)^/(2(i„-i„_i)) 

/ {{tn-tn-l)\l'2^T^{tn-tn-l))dXn 

is nonzero for all an,f3n, since < c„-i. Hence, Dini's theorem holds 

and we get that (3n = Pn{ctn) is a continuous and differentiable function. If 
Sn = On + ^1 the error on a„ is propagated on /?„ as 

where 

A/3„ = /?„(«„ + 5) -/3„(a„)~0(5). □ 

Theorem 4.3. The error of the PLMC method at the discretization 
knots tn, n = 1, 2, . . . , is 

|ePLMC| ^ |^(^^) _ c(t„)| ~0(max(5,/i2)) 

w/ien ai is affected by an error of the order of 6 >0. 

Proof. When ai is affected by an error of the order of 5, by Lemmas 4.1 
and 4.2, we get that 

\b{ti) - c{ti)\ < - ci(ti)| + |ci(ti) - c(ti)| 

^0{h^) + 0{6) 

~0(max(5, /i^)). 

By induction, let — c(t„_i)| ~ 0(max((5, /i^)), then 

= \b{tn-l) - Cn-1 + (/3n " Pn)h\ 
< \b{tn-l)-Cn-l\ + \{Pn-Pn)h\ 

~ 0(max((5, /i^)) + 0(max((5/i, /i^)) 
~0(max(5, h^)) 

and we obtain 

\b{tn) - C{tn)\ < \bitn) - C„(t„)| + |c„(t„) - c(t„)| 

~0(/i2) + 0(max(5, /i^)) 

~0(max(5,/i2)). □ 

Remark 4.1. If the boundary b{t) is not monotone nor convex (or con- 
cave) but is sufficiently well behaved, one can proceed with a similar reason- 
ing on each monotonicity and convexity interval. Hence, choosing km a suit- 
able way, the stepwise boundary has still an error jej^^'^^l '-^ 0{max{S,h'^)). 
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4.2. Error due to the Monte Carlo approximation. Disregarding the error 
related to the numerical integration of the l.h.s. of (3.1) at the first step of 
the algorithm that can be well controlled with a careful use of numerical 
approximations, we focus on the next steps when the multiple integrals of 
the l.h.s. of (3.2), 



rtn+i 


rc„ 


■f ( 




' — oo 


J — oo \ 



F{Pn+l)= / •••/ [Mt\tn,Xn 



i,Xj_i) dxi ■ ■ ■ dxndt, 



X 

1=1 



are evaluated via Monte Carlo method. Using the Law of Large Numbers, 
we approximate the expectation on the l.h.s. of (3.5) with its sample mean. 
Hence, for a fixed confidence level a, we get 



P 



1 



n 

X 

=1 



i^(A 



< 6n \ > a. 



where Xnj ~ N{0, ti — tj-i), i = 1, . . . , n, j = 1, . . . , M, are independent ran- 
dom variables. Letting i = 1, . . . , n, j = 1, . . . , M, be a sample of X^^j, 
with an accuracy 5^ in the computation of the integral, we obtain a confi- 
dence interval for (3n+i at the same level a: 

( \ 

F''[j^T.Hi(3n+i;x^,,) 

n N 

X n^(^ij ^ Ci)e-2('='-^^>^)('='-i-^^-i.^)/(*»-*^-i) ±5« 

1=1 / 

Remark 4.2. The computations necessary to get the Monte Carlo eval- 
uations are not excessively expensive. Hence, we can choose a very large 
value for the size M for the number of simulations involved in the integral 
estimation. This allows to make this error negligible with respect to the 
error determined by the piecewise linear approximation of the boundary, 
but as a consequence of the use of the Monte Carlo method, our results are 
characterized by a reliability a. 
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5. An approximation by the nonlinear Volterra integral equation. The 

algorithm discussed in the previous section is reliable and easily imple- 
mented, but it is computationally expensive since the Monte Carlo method 
requires long times of computation. In this section we therefore consider 
an alternative approach of pure numerical nature which is computationally 
more attractive. Since this method is based on the approximation of a non 
linear Volterra integral equation, it will be referred to as VIE (Volterra In- 
tegral Equation) method. 

Let us consider the integral equation [cf. Peskir (2002)]: 



(5.1) ^ 



where ^{x) = 1 — $(x) and <I>(x) = viz) dz is the standard normal dis- 
tribution. Equation (5.1) is a Volterra integral equation of the first kind in 
/;,, but it is a nonlinear Volterra integral equation of the second kind in h 
and its kernel is nonsingular in the sense that it is bounded. Moreover the 
nonlinear functions involved in the equation are bounded and ^ is invertible. 
These features allow the development of the following numerical algorithm 
that approximate (5.1) in a simple and reliable way. 

We numerically solve this equation to evaluate the approximate value h* 
of h at the knots = ih for i = 1, . . . , n, where h = t/n (and t > is given 
and fixed). To this aim we follow the original idea of Volterra [see, e.g., Linz 
(1985), Chapter VII], that is, we approximate the integral on the l.h.s. of 
(5.1) with the Euler method 

getting a nonlinear system of n equations in n unknowns b*{ti), . . . ,6*(t„). 
The ith equation of (5.2) for i > 2 makes use of the values b*(ti),. . . , 6*(tj_i) 
found in the preceding steps. Hence, equations (5.2) can be solved itera- 
tively using the iterative middle point method [cf. Atkinson (1989)], which 
then gives approximate values for the unknown boundary b at the points 

^1 ) • • • 1 • 

We recall that the local consistency error for (5.1) for a generic discretiza- 
tion method of the integral is [cf. Linz (1985)] 

(5.3) 



j=0 



3)^ 
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where coij are the integration weights of the adopted discretization schema. 
In (5.2) we use the Euler method for which we have tOij = 1 and Wjo = 0, for 
each j = 1, . . . , i, i = 1, . . . ,n. Since maxo<i<n 6{h, ti) = 0(h), the method is 
consistent of order 1. 

Remark 5.1. Note that the system (5.2) is triangular and this makes it 
especially efficient if one wishes to compute b* at the next knot, when it is 
known in the previous ones. Moreover, unlike the PLMC method, here the 
knowledge of 6(0) is not required since we use the "forward" Euler method, 
also known as right-hand rectangular rule [cf. Atkinson (1989)]. 

Remark 5.2. In spite of the Euler method, we could use the extended 
trapezoidal formula [cf. Abramowitz and Stegun (1964), page 885, formula 
25.4.2] with weights uJio = cou = 1/2 for each i = 1, . . . ,n and tUij = 1 for each 
J = 1, . . . , i — 1, i = 1, . . . , n, to approximate the integral in (5.1). In this case 
(5.2) is replaced by 

.|.(M^).(,,. 

+ ^fb{ti)h (i = l,...,n). 

This approximation is consistent of order 2, but in the general case re- 
quests the knowledge of fe(0). However, when /^(O) = the first term in the 
r.h.s. of (5.4) vanishes and we obtain again a triangular system independent 
from the knowledge of 6(0). 

Remark 5.3. The unnecessity of the knowledge of 6(0) is a numerical 
advantage, but it hides some problems connected with the initial point ti. 
To illustrate this, note that the algorithm uses the approximation (5.2) that 
for i = l reads 

(5.5, Kl^)4«"'"' 

using that ^(0) = 1/2 (when 6 is smooth). Recall that ti = t/n so that 
ti — > as n — > oo. Taking then b{t) = \j2t log(l/t), for example, it is easily 
verified that lim^^oo 2^'(6(ti)/v%)/ti = 0, while the result of Peskir (2002) 
implies that /;,(0+) = +oo. Thus, the approximation (5.5) in this case fails 
in a neighborhood of zero. A closer look shows that it is valid if /fe(0+) = 
(under certain mild regularity conditions). It follows that the more /fe(0+) is 
away from zero, the less accurate (5.5) becomes in a neighborhood of zero. 
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Remark 5.4. To improve the boundary estimation for small times, we 
can consider the flux equation introduced in Ricciardi, Sacerdote and Sato 
(1984) and its approximation as t — > to get 



(5.6) m^^^e^p 



bjtf 
2t 



t<e, 



where e is small enough to make it acceptable to approximate the r.h.s. of 
the flux equation with its first term. It is an implicit nonlinear equation 
in terms of b{t) that can be numerically solved. Since this equation loses 
reliability as the time t increases, it is necessary to control the validity of 
this equation for the value of t of interest. Substituting the solution b{t) of 
(5.6) in the complete equation (2.8) in Ricciardi, Sacerdote and Sato (1984), 
we obtain an estimator fb{t) of the FPT density function. The evaluation 
of the relative error between this estimator and the known value of the 
density function allows to guess the time interval where this approximation 
is reliable. 



6. VIE method error estimation. Nonlinear integral equations are con- 
sidered in Linz (1985), but the integral equation (5.1) differs from nonlinear 
integral equations studied there due to the expression on the l.h.s. Indeed, 
it contains the unknown function in an implicit way through a further non- 
linear function. However, working in analogy with the methods used in Linz 
(1985), we can prove the convergence of the VIE method of Section 5. To 
this end, we recall a convergence theorem used by Linz (1985). 

Theorem 6.1. Let the sequence S^o, ^i,... satisfy 

n-l 

(6.1) \U<^T.\^i\+Bn, n = r,r + l,..., 

i=0 

where ^4 > 0, |-B„| < B and assume that it exists a constant r] > such that 

r-l 

(6.2) El^^l^^- 

i=0 

Then 

(6.3) \^n\<{l + Ar-'-{B + Ar^), n = r,r + l,.... 

We now can prove the convergence of the proposed method. 

Theorem 6.2. The error e^^^ of the VIE method at the discretization 
knots tn, n = 1, 2, . . . , is 

(6.4) \el'''\ = \h{tn)-h*{t^)\^0{h) 
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if the integral in (5.1) is approximated via the Euler formula. 

If fb{0) = and the extended trapezoidal formula is used for the integral 
in (5.1 ), the error becomes 

(6.5) \el'''\ = \b{tn)-h*{tr,)\^0{h^). 

Proof. We split the proof of (6.4) in two parts, the first one concerning 
the first discretization interval and the second one concerning a generic step 
n. 



Step 1. Choosing t = ti and subtracting (5.2) from (5.1), we get 



(6.6) 



b*iti 



Using the inverse function ^' ^, we obtain 



(6.7) 



b{ti 



-1 



b*iti 



+ S{h,ti 



Substituting the Taylor expansion of ^' ^ around y = ^( ^J^^ ) in this last 
equation, we get 



(6.^ 



dy 

+ V2^6{h,ti)e^^/^ 



y=^{b*{ti)/y^)+e5{h,ti) 



where 9 G (0, 1) [cf. Abramowitz and Stegun (1964)] and 



(6.9) 



b*{ti) 



+ e5{hM) 



with L < oo. 

To prove (6.9), note that if 5{h,ti) > 0, one has 



(6.10) 



b*iti 



+ 9S{h,ti) > ^ 



b*{ti) 



and, due to the decreasing monotonicity of ^ ^, one gets 

b*{ti) 



(6.11) 



<Li. 



Elsewhere, when 6{h,ti) < one has 

'b*{t,y 



20 

hence, 
(6.13) 



C. ZUCCA AND L. SACERDOTE 



b{t 



Choosing L = max{Li, L2), we get (6.9). Hence from (6.8) and (6.9) we get 



^VIEi 

^1 I 



\b{h)-b*{ti)\ < \V2^S{h,ti)e^^/^\<k\6ih,ti)\=Oih), 



since the Euler method is consistent of order 1. 



Step n. Choosing t = tn and subtracting (5.2) from (5.1), we get 
bitnV 



(6.14) 
where 

(6.15) 



n = 2, 3, . 



7n = lnih,ti, ...,tn) 

+ 5{h,tn), n = 2,3,.... 



b*{tn)-b*{tj) 



Mimicking the procedure used for case n = 1 , we apply the inverse function 
^'^^ to (6.14) and we expand in Taylor series the resulting r.h.s. around 



^ = ^-(^1111). Thus, we g( 



;et 



b{tn) b*{tn) , d^-^y) 



(6.16) 



b*{tn) 



dy 

+ \/27r7„e 



s/=i'(fe*{t„)/v^)+e7„ 



where 9 G (0, 1) [cf. Abramowitz and Stegun (1964)] and 



(6.17) 



with Mn < 00 for each n = 1, 2, 

The last inequality can be proved analogously to (6.11) and (6.13) in 
Step 1. Hence, from (6.16) we get 



(6.18) 



,VIE| 



\b{tn) - b*{tn)\ < I V2^7ne*^"/'| < k\jn\. 



Let us now split the term 7„ as the sum of a contribution of the ac- 
cumulated error due to the previous steps and a contribution due to the 
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consistency of the numerical approximation of the integral in (5.1). To this 
goal, we bound the module of 7„ using the triangular inequality to get 



(6.19) 



hn\<hY,fb{tj)^ 

+ \6{h,tn)\ 



n-1 
+ \5ih,tn)\, 



bjtn) - bjtj 



b*{tn)-b*{t,) 



where the last equality uses the definition of ^ and the fact that the last 
term in the sum is zero. Bounding the integrand in (6.19) with 1, we get 



n-l 



(6.20) |7n|</iE 



^'^''^ (kn-k7^^l) + |5(/^,tn)| 



1 W 27r(t„ — t 



and (6.18) becomes 



n-l 



^/2TTtnMh 



■J) 



VIE: 



(6.21) 



j=i Y 27r(t„ — tj] 

n-l 



+ V^-KtnMh 



fbitj 



]=i Y 27r(t„ — tj ^ 
where M = max„{exp(M„)} and, hence. 



\eJ^^\ + V2^M\5{h,tn) 



n-l 



fbitj) 



(6.22) 
If 

(6.23) 
we get 
(6.24) 



j=i h , 



n-l 



V*n tj 



h< 



.-i<^Ei.rK^i.(M.)i, 



1 - hk 
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where F = sup,- „ —J=^l== . 

Applying the convergence Theorem 6.1 with 

MhF 

(6.26) B„ = ^^^\S(h,t„)\<B 

1 — hk 

and noting that by inductive argument one has 

n-l 

(6.27) T.\^T\<V. 
we finahy obtain 

(6.28) \el'^^\<(^^^\5{h,tn)\2T^M + MhFY,^\ef^ 

Recalhng that the method is consistent of order 1, part 1 of the theorem 
is proved. Using equation (6.21) in the case of the extended trapezoidal 
scheme, one easily proves formula (6.5). □ 

Remark 6.1. Note that the error of the method is dominated by the 
consistency error. 

Remark 6.2. To prove Theorem 6.2, we simply use the monotonicity 
properties of the function ^ . Hence, the method can be easily extended to 
different diffusion processes, simply substituting ^ in (5.1) with the survival 
function of the considered process. 

7. Examples. In this section we check the stability of the algorithms 
presented in Sections 3 and 5 by means of some examples where a closed 
form solution is available. We also show other examples where the solution is 
numerically evaluated. First we apply the algorithms of Sections 3 and 5 for 
a Daniels boundary h for two sets of parameters (Section 2). Later we apply 
them to the case of an oscillating boundary with different parameters. In 
these later cases the FPT density function has been numerically estimated 
via the Buonocore, Nobile and Ricciardi algorithm (1987). We consider the 
mean square deviation as an index of the goodness of the two methods 

1 " 

(7.1) 4i) = ^(6(t.)-6«(t^.))', i = l,2, 

where b^^^ denotes the approximating boundary determined by the PLMC 
algorithm of Section 3 and 6^^^ denotes the approximating boundary deter- 
mined by the VIE algorithm of Section 5. 

We apply the two algorithms to the following two cases: 
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1. Daniels boundary with parameters a = = 0.5,7 = 0-5 and a = 
1,(3 = 1,7 = 0.5. We first apply the PLMC algorithm in the interval [0,2] 
with a discretization step h = 0.2 and we compute the integrals by a Monte 
Carlo method using 10^ simulations. Under these conditions we estimate 
the mean square deviation, obtaining Un ^ = 9.4 • 10~^ and an^ = 3.4 • 10~^ 
respectively. 

We also apply to these examples the VIE algorithm of Section 5 with 
discretization step h = 0.01 and the resulting mean square deviations are 
a'k^ = 4.3 • 10-5 and a^n^ = 4.6 • 10^^ respectively. 

2. Oscillating boundary 

(7.2) b{t) = a + (3 cos{-ft) 

with parameters a = l,/3 = 0.5,7 = 2 and a = l,/3 = l,7 = 2 respectively. 
Repeating the PLMC algorithm under exactly the same conditions as above, 
we get the following mean square deviations: = 6.6 • 10"^ and = 
6.4 • 10"^ respectively. The application of the VIE algorithm to the two sets 

(2) (2) 

of parameters gives dn = 4.9 • 10 and an = 3.4 -10 . 

These results confirm the reliability of the methods. The four cases are 
illustrated in Figures 2, 3, 5 and 6 where the exact boundary shape (left) is 
compared with the approximating ones obtained by means of the algorithms 
of Section 3 (center) and Section 5 (right). 

Remark 7.1. With reference to the PLMC method, we underline that 
the goodness of the approximation does not depend only on the discretiza- 
tion step h, but also on the probability mass for the boundary to be crossed 
in the discretization subinterval of length h. It is thus recommendable to 
avoid excessively small discretization steps. However, as shown in our exam- 
ples, quite large values of h give small value for the mean square deviation. 

Remark 7.2. A comparison between the exact and the approximating 
boundary obtained from the PLMC algorithm in Figures 2 and 3 with the 
first-passage density function in Figure 4 shows the rise of larger oscillations 
as t increases. This fact can be explained by observing that as t increases 
the probability of crossing the boundary on each discretized interval becomes 
smaller and the Monte Carlo method is subject to a larger relative error. A 
further improvement of the method could be to use an adaptive step built 
on a constant probability mass of the first-passage density function on each 
interval. 

Remark 7.3. In Figures 2 and 3 we observe that the approximating 
boundary obtained from the VIE algorithm of Section 5 has a large error 
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0.5 1 1.5 2 0.5 1 1.5 2 1 2 

\ t t 

Fig. 2. Daniels boundary with parameters a = = 0.5,7 = 0-5 fl^ft) compared with the 
approximating ones obtained by means of the PLMC method (center) and the VIE method 
(right). 



for short times. This result can be explained by the rough approximation 
in the first steps of the present method and is related to the use of the 
Euler method. In Remark 5.4 we introduced an improvement of the method 
for small times. In Figure 7 we compare the approximation obtained via 
the VIE method (dash dot line) and two different corrections via (5.6) in 
a special case of Daniels boundary (full line). The stars indicate the values 
of the boundary solution of (5.6), while the two other curves substitute the 
first one or two values of h{t) with the solution of (5.6). We observe that the 
new boundary estimates are more reliable. 

Remark 7.4. The evaluation of the difference between the exact bound- 
ary and the approximating one, obtained with the PLMC method, shows 
that the approximating boundary oscillates around the real one. 




Fig. 3. Daniels boundary with parameters a = 1,13 = 1,7 = 0.5 (left) compared with the 
approximating ones obtained by means of the PLMC method (center) and the VIE method 
(right). 
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0.5 1 1.5 2 0.5 1 1.5 2 

t t 

Fig. 4. The first-passage density function through a Daniels boundary with parameters 
a = 1, /3 = 0.5, 7 = 0.5 (left) and with parameters a — 1, (3 = 1,^ = 0.5 (right). 

8. The exponential case. Let us apply our algorithms to numerically 
solve the original Shiryaev's problem, that is, to approximate the boundary 
corresponding to an exponential first-passage density function. 

Let us thus consider a first-passage density function fh{t) = e~* for t > 
corresponding to the exponential distribution with parameter A = 1. In 
order to apply the PLMC method, we note that this particular choice of 
distribution requires to study a first-passage time T), such that fbiO+) > 0. 
As indicated in Section 2.2, this implies that the boundary b should be an 
upper function for W that vanishes at zero so that b'{0+) = +oo. 

Making use of (2.14) and introducing more generally the following nota- 
tion: 

(8.1) « = /(,(0+)>0, 




t t t 



Fig. 5. Oscillating boundary with parameters a — l,l3 = 0.5,7 = 2 (left) compared with 
the approximating ones obtained by means of the PLMC method (center) and the VIE 
method (right). 



26 



C. ZUCCA AND L. SACERDOTE 



to 




Fig. 6. Oscillating boundary with parameters a=l,/3 = l,7 = 2 (left) compared with the 
approximating ones obtained by means of the PLMC method (center) and the VIE method 
(right). 



we obtain an approximation of the boundary b in the neighborhood of zero 
by choosing 

(8.2) c = -21og(\/4^K). 

On the first interval we approximate the boundary b by (2.13), where c is 
given by (8.2). On the successive intervals we simply apply the algorithm 
with 



i.3) 



C2 = gih) = \2tl log(l/tl) + h loglog(l/tl) + Cti. 




0.55 



0.45 



0,02 0.04 0.06 O.OS 0.1 0.12 0.14 0.16 0.1 S 
t 

Fig. 7. Daniels boundary with parameters a — 1, (5 = 0.5,7 = 0-5 (full line) compared 
with the approximating ones obtained by means of the VIE method (dash dot line) and 
by means of the VIE method modified for small times (dashed lines) substituting the first 
one or two values of b(t) with the solution of (5.6). The stars indicate the values of the 
boundary obtained as a solution of (5.6). 
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0,5 1 1.5 2 0.5 1 1.5 2 

t t 



Fig. 8. Numerical evaluations of the approximating boundaries for an exponential first- 
passage density function with parameter A = 1 obtained by m,eans of the PLMC method 
(left) and the VIE method (right) with a discretization step h = 0.01. 



In Figure 8 we plot the boundaries corresponding to an exponential first- 
passage density function with A = 1 obtained by the two different algorithms 
described in the previous sections; the plot on the r.h.s. is obtained applying 
the PLMC method with the step h = 0.01 and with a number of simulations 
for the Monte Carlo method equal to 10^, while the plot on the l.h.s. is 
obtained applying the VIE method with discretization step h = 0.01. 

A heuristic confirmation of the stability of the algorithms is given by the 
nearness of the two curves. Furthermore, an intuitive reasoning about the 
shape of the boundary confirms our result. To obtain an exponential first- 
passage density function, a large part of the sample paths should cross the 
boundary very soon but, as the time goes on, other sample paths, that have 
not yet reached the boundary, should be able to cross it too. This makes 
intuitive the fact that the boundary must decrease in order to be reachable 
by the sample paths whose initial trend was negative. 
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